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ABSTRACT 

The  development  of  microelectronics  has  brought  into  being 
large  capacity  digital  memories  in  a  small  package.    In  the 
foreseeable  future  even  more  advances  can  be  seen  in  this 
trend.    Therefore  the  use  of  digital  computers  in  control  sys- 
tems will  play  an  even  larger  role  than  today. 

This  work  involves  a  fourth  order  system  to  simulate  the 
control  and  dynamics  of  a  missile.    Proportional  navigation 
is  used  as  the  guidance  method.    Studied  are  the  effects  of 
applying  different  controls  which  are  considered  best  from  a 
computer  study  and  the  effects  of  applying  digital  filtering 
methods.    Although  these  studies  were  applied  to  a  specific 
problem,  an  attempt  is  made  to  keep  the  discussion  general 
in  order  that  the  methods  may  be  considered  for  other  problems. 
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CHAPTER  I 
INTRODUCTION 

If  it  is  desired  to  hit  a  target  with  a  missile,  some  form  of 
guidance  must  be  employed.    This  guidance  would  be  selected 
weighing  the  desired  objectives  and  an  estimate  of  the  target 
capabilities  against  the  amount  of  money  and  missile  space 
available.    As  suggested  in  the  introduction,  the  advances  in 
computer  size  and  capabilities  is ,  and  will  in  the  future ,  con- 
tinue to  change  the  weighting  of  these  factors .    Thus  in  the 
foreseeable  future  it  may  be  possible  to  design  a  model  of  the 
system  which  will  be  in  sufficient  detail  to  predict  the  respon- 
se of  the  real  system  to  a  high  degree  of  accuracy.    Using  the 
micro-computers,  this  model  could  then  become  a  part  of  the 
control  system  in  the  missile.    Such  a  control  system  might  be 
described  in  the  block  diagram  in  Figure  1-1 . 

In  Figure  1-1  the  missile  measures  and  reacts  to  a  signal 
which  depends  on  the  type  of  guidance  employed.    Such  a 
signal  might  depend  on  the  rate  of  change  in  the  missile's 
radar  seeker  head  which  is  tracking  the  line  of  sight  from  the 
missile  to  target.    The  missile  will  react  to  this  signal, pro- 
ducing through  its  dynamics  a  change  in  the  direction  of  its 
velocity  vector,  or  the  missile  will  take  some  other  action 
called  for  by  the  particular  guidance  law  being  used.    The 
model  receives  the  measured  signals  and  considers  the 
measurable  states  from  the  missile  dynamics.    Since  all  the 
states  of  the  model  are  measurable  and  the  digital  model  in 
itself  generates  relatively  little  noise,  the  model  becomes  a 
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Figure  1-1.  Block  Diagram  of  a  Dynamic  System, 
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state  predictor  or  estimator. 

Since  all  the  states  of  the  model  are  available,  a  control 
law  considering  all  the  states  can  now  be  designed.    This 
control  law  is  applied  in  the  controller  to  guide  the  missile. 

If  the  desired  objectives  of  the  system  can  justify  the 
above  approach,  the  system  performance  of  the  modified  system 
should  be  an  improvement  over  most  conventional  guidance 
systems . 
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CHAPTER  II 

PROPORTIONAL  NAVIGATION 

Adler  in  [9)]  defined  proportional  navigation  as, 

A  course  in  which  the  rate  of  change  of  missile  heading  is 
directly  proportional  to  the  rate  of  rotation  of  the  line-of- 
sight  from  the  missile  to  the  target. 

y=Kcr  (1) 

where 

yis  the  angular  rate  of  change  of  the  missile  velocity  vector. 

(jis  the  angular  rate  of  change  of  the  line-of- sight. 

K  is  a  constant,  typically  between  three  and  five. 

To  gain  a  better  understanding  of  proportional  navigation, 
some  of  the  simpler  forms  of  line-of-sight  guidance  will  be 
examined.    Pursuit  course  (sometimes  called  pure  pursuit)  always 
aims  the  missile  directly  at  the  target  along  the  line-of-sight. 
This  method  will  always  end  in  a  tail  chase,  even  if  the  missile 
is  launched  head-on  with  the  target.    High  missile  acceler- 
ations are  required.    It  can  be  shown  from  the  equations  of 
motion  that  if  the  missile  to  target  speed  ration  exceeds  two, 
the  final  missile  turning  rate  will  approach  infinity.    Thus 
pursuit  course  guidance  may  not  be  very  satisfactory  although 
it  is  very  simple  to  implement. 

The  next  step  is  to  lead  the  line-of-sight  by  an  angle 
which  is  a  function  of  the  target  velocity.    This  type  is  called 
constant  bearing  course,  and  is  achieved  by  aligning  the 
relative  missile  to  target  velocity  vector  with  the  line-of- 
sight.    That  is  to  say  the  line-of-sight  maintains  a  constant 
direction  in  space.    It  can  be  shown  that  this  method  can 


14 


cause  a  missile  to  follow  a  target  to  a  collision  even  if  the 
target  is  maneuvering.    This  method  however  requires  an  in- 
stantaneous correction  to  each  change  in  the  line-of-sight. 

Consider  the  constant  bearing  trajectory  shown  in 
Figure  2-1.    In  the  figure  the  line-of-sight  is  shown  for 
several  samples  made  by  the  radar.    As  long  as  the  missile 
remains  on  the  collision  course,  the  angular  rate  of  change  of 
the  line-of-sight  will  remain  zero  and  no  control  action  is 
required.    Any  rotation  of  the  line-of-sight,  indicating  a  de- 
parture from  the  collision  course  will  be  detected  by  the 
missile's  radar.    The  missile  using  proportional  navigation 
will  turn  at  a  rate  proportional  to  this  rotation  in  a  direction 
to  reduce  the  line-of-sight  rate  and  return  to  a  constant 
bearing  collision  course. 

In  the  Figure  2-2  above,  the  lead  angle  /Sis  defined: 

j3=<7-y  (2) 

This  angle  (/3)  is  independent  of  any  reference  used  to  define 
o  and  y.    If,  as  in  the  discussion  above,  we  could  achieve  a 
constant  bearing  course,  /??  would  remain  a  constant,  and 
therefore  ^ would  remain  zero. 
Thus  for  a  constant  bearing  course , 

£=a-y=0  (3) 

Since  the  missile  cannot  react  instantaneously,  /Jis  not 
always  zero.    Thus  for  any  instant  of  time  (1)  may  be  written 
as , 

<i=hy  (4) 

If  /Sis  equal  to  zero,  the  missile  is  on  a  constant  bearing 
collision  course.    Therefore  the  object  of  the  control  system 
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Figure  2-1.  Constant  Bearing  Trajectory. 


Figure  2-2.  Proportional  Navigation  Geometry, 
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Figure  2-3.   Missile  Guidance  Inter  Loop. 
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in  proportional  navigation  is  to  minimize  £and  if  CT becomes  a 
constant  or  zero  to  drive  j8to  zero. 
Guidance  and  Control 

A  typical  homing  missile  inter  guidance  loop  is  shown  in 
Figure  2-3.    The  following  discussion  will  be  limited  to  only 
coplanar  motion  of  the  missile  and  target.    A  constant  missile 
and  target  speed  is  also  assumed  throughout  the  rest  of  this 
paper. 

In  Figure  2-3,  the  radar  tracks  the  target  using  a  servo 
loop  to  keep  the  antenna  on  the  line-of-sight.    A  delayed  sig- 
nal which  is  a  measure  of  cr,  is  obtained  from  the  radar  system 
A  voltage  which  is  proportional  to  the  rate  of  change  of  this 
signal  is  obtained.    This  voltage  is  compared  with  the  output 
of  the  dynamic  system  y,  and  the  error  signal  is  obtained. 
This  error  signal  is  sent  either  to  the  hydraulic  system  to 
rotate  the  control  surfaces  or  to  a  hydraulic/valve  system 
which  controls  side  thrust  jets,  depending  on  the  type  of 
missile  being  used. 
Kinematics 

As  might  be  suggested  by  the  name  inter  loop,  there  is 
also  an  outer  loop  as  shown  in  Figure  2-4.    In  the  Figure, 
the  guidance  and  control  block  is  shown  as  in  Figure  2-3. 
The  effect  of  yon  the  line-of-sight  and  the  effect  of  the  tar- 
get motion  on  the  line-of-sight  are  represented  in  the  missile- 
target  kinematics  block. 

This  outer  loop  is  not  as  accessible  to  measurement 
and  control,  but  the  understanding  of  its  effects  is  necessary 
in  order  to  perform  a  simulation  of  the  entire  system.    From 
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the  geometry  of  Figure  2-2,  the  following  equations  may  be 

written: 

R  =  Vt    COS(a)  -  Vm  COS(/3)  (5) 

Rcr=  -Vt  SIN(ct)  +  Vm  SIN(0)  (6) 

where 

Vt    =  Target  speed  (assumed  Constant) 

V     =  Missile  speed  (assumed  Constant) 
m 

R  =  the  line-of- sight  distance 

R  =  the  rate  of  change  of  the  line-of-sight. 
Then  equations  (5)  and  (6)  may  be  rewritten: 

R  =  Vt  COS(ct)  -  Vm  COS(<7  -  y)  (7) 

Ra=  -Vt  SIN(cr)  +  Vm  SIN(a-  y)  (8) 

Differentiating  (8)  with  respect  to  time: 

Ra+Ro=  -VtCOS(a)  <r+  Vm  COS(CT-  y)(o-  y) 

=  -Ra-  VmCOS(a-y)y  (9) 

Therefore 

2Ra+  Ra=  -V     COS(a  -  y)  y  (10) 


• 

R  = 

-R/T 

and 

Vr 

=tei 

where  T  is  the  time  to  go  to  the  target.    Then  equation  (10) 
may  be  represented  by: 

(-2R/T  +  RS)  <J  =  -Vm  COS(a  -  y)  y  (11) 

thus 

°°/y=-Vm  COS(o-y)  (12) 


R/T  (TS  -  2) 
COS 
R/T 


Where  V   K  COS(a-  y)  is  called  the  effective  navigation  constant. 

R/f 
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Lil 


Since  Vm  is  constant  and  variations  in  GOS(a-  y)  and  Vr  are 
small,  this  term  is  assumed  constant  and  values  of  four  to 
six  are  considered  best. 
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CHAPTER  III 
FILTER  AND  CONTROL 
Plant  Description 

A  linear,  time-invariant  dynamic  system  is  described  in 
the  flow  graph,  Figure  3-1 . 
The  transfer  function  of  Figure  3-1  may    be  stated: 

cout 


.n-1 


X_,+  (s)         \s11-1  + +  b2s  +  b1 


X, 


;n+a„sI1-i  + 


(1) 


(s)  s11  +  ans  "  -1  +  .   .  +  dL2s  +  a1 

Equation  (1)  may  be  written  in  one  differential  equation  using 
matrix  notation. 

X=FX+DU,C  =  BX  (2) 

Where 

X  is  the  state  vector  (n  x  1) 

U  is  a  scalar  of  the  system  inputs  and  controls 

F  is  a  matrix  of  constants.    For  the  system  described  in 
Figure  3-1,  F  is  defined  in  equation  (3). 


F  = 


0 

1 

0 

\                   • 

0 

0 

0 

1 

• 

0 

• 

0 

• 

0 

• 

0 

*               • 

1 

• 

0 

0 

0 

0 

0 

1 

-a 
L    i 

~*2 

• 

»                   • 

-a 

(3) 


D  is  a  matrix  of  constants  (n  x  1)  described  in  equation 
(4)  for  the  system  in  Figure  3-1 . 
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Figure  3-1.   Flow  Diagram  of  the  system. 
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D  = 


(4) 


B  is  an  (1  x  n)  row  vector,  described  in  equation  (5)  for 
the  system  described  in  Figure  3-1 . 

B=  [bj  b2  .    .    .   bn]  (5) 

The  solution  to  (2)  is  given  by: 

X(t)  =$(t  -  t  )  X(tQ)  +     Jt^  (t  -  tx)  D  dtx  U(t0)      l        (6) 


-1 


Where 

*(t  -  tQ)  =    <=C_1(SI  -  F) 

Solving  the  differential  equation  (2)  for  the  discrete  solution, 
X(k  +  1)  =  $X(k)  +  DEL  U(k)  (7) 

Where 


<£ 


=  eF(t-  tQ) 


DEL 


=    /U^'-^D 


to 


dt1 


(8) 
(9) 


Note  that  k(t  -  t  )  is  represented  by  k. 

Titus,  in  reference  [10]  developed  a  digital  computer 
program  called  PHIDEL.    This  program  provides  a  solution  to 
equations  (8)  and  (9)  and  will  be  used  as  a  subroutine  to 
obtain  <i>and  DEL.    A  listing  of  the  PHIDEL  program  is  provided 
in  Appendix  one. 
Plant  Control 

It  is  desired  to  design  a  set  of  feedback  coefficients  to 
control  the  plant  in  accordance  with  a  selected  performance 

1  U(t)  is  held  constant  over  the  interval  (t  -  tQ)  and  is 
equal  to  U(tQ). 
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index.    This  set  of  coefficients  will  modify  the  plant  states  and 
in  effect  will  cause  a  movement  of  the  Z-plane  poles.    This 
effect  must  be  weighed  against  the  feedback  gains  derived  for 
the  desired  performance  index. 

Let  J(N)  be  defined  as  the  performance  index  of  a  discrete 
sample-data  system  which  will  be  minimized  with  respect  to 
U(K),  the  control  for  the  plant. 

N  *  9 

J(N)  =  Min    E     [X^K)  Q  X(K)  +rUz(K-l)]  (10) 

U(K)   K  =  M 

Q  is  a  (n  x  n)  positive-definite  symmetric  constant  matrix, 

and  R  is  a  positive  scalar  constant. 

The  principle  of  optimality  states  that  any  portion  of  an 

optimal  trajectory  is  also  an  optimal  trajectory.    Therefore 

equation  (1)  may  be  rewritten: 

J(N)  =  Min      [X^N)  Q  X(N)  +  r  U2(N  -  1)  +  J(N  -  1)]        (11) 
U(N) 

If  the  gradient  of  J(N)  is  taken  with  respect  to  U(N  -  1)  and 

set  equal  to  zero,  then  an  optimal  control,  U°(N  -  1)  may  be 

determined.    It  is  noted  that  J(N  -  1)  is  not  a  function  of 

U(N  -  1). 

Using  the  above  arguments,  Titus  developed  the  following 

algorithms: 

AMk  +  1)  =  -&  P(k)  3- 

AlP(k)  A+  r  (12) 

P(k)  =  <&  (k)P(k  -  1)  $(k)  +  Q  +  rA(k)Al(k)  (13$ 

Mk)  =  $+  AAt(k)  (14) 

Using  the  recursive  relations  of  (12),   (13)  and  (14),  Titus    [10] 

developed  program  OPTCON .    This  method  is  also  discussed 

in  references  [3]     and     [4]  by  Titus,  Strum  and  Demetry,  and 
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in  reference  [6  ]by  Ogden.    A  fortran  listing  of  the  program  may 
be  found  in  Appendix  one . 
Digital  Filter 

In  the  control  discussion,  the  assumption  was  made  that 
all  the  states  of  the  dynamic  system  are  observable  states.    This 
condition  is  certainly  not  always  true  and  often  the  observable 
states  may  be  measured  only  at  the  cost  of  some  ambiguity  due 
to  the  measurement  noise.    The  inputs  to  the  system  such  as 
the  radar  seeker  discussed  in  Chapter  Two  will  often  be  re- 
ceived in  a  noisy  environment. 

Kalman  in  references  [l]   and    [2  ]  discussed  these  prob- 
lems and  presented  the  theory  for  the  desired  filter.    Schmidt 
(5  ]  ;  Titus,  Demetry  and  Strum    [4]  and  Jardine  [7]   have  dis- 
cussed this  problem  and  developed  methods  of  implementing  this 
problem  on  the  digital  computer.    In  an  effort  to  maintain  claftty, 
some  of  these  developments  will  be  presented. 

The  digital  filter  will  provide  a  best  estimate  of  all  the 
states  by  weighing  the  past  information  with  the  present  obser- 
vable states.    This  weighting  is  performed  with  the  knowledge 
of  the  environment  (excitation)  and  the  measurement  noise . 
Although  in  many  cases  the  noise  is  very  difficult  to  describe, 
the  assumption  that  white  noise  is  present  is  in  general  true. 
Thus  if  white  noise  can  be  discriminated,  the  integrity  of  the 
measured  signals  will  be  improved. 
Filter  Design 

For  a  single  variable,  the  desnity  function  f(x)  for  a 
gaussian  distribution  is  given  by: 

f(x)  =  1  e-&x  r-  u)2 

(2ir)i      o  —^2—  dS) 
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where 

£g(x)f(x)dx  =  E[g(x)]  (16) 

where  E  implies  expected  value, 
and 

5*f  (x)dx  =  1  (17) 

Assume  that  the  systems  measurement  and  excitation 
noise  are  white  noise,  where  white  noise  has  a  gaussian  dis- 
tribution and  is  spread  uniformly  over  all  frequencies. 

Noting  that  the  input,  U(k)  to  the  plant  is  independent  of 
U(k  -  1),  X(k)  is  independent  of  X(k  -  1)  and  using  the  joint 
probability  property  of  random  independent  samples,  a  density 
function  for  the  state  system  may  be  written: 

F(X)=  1  e-^X-X^p-^X-X]  (18) 

(2  7r)Vz  pi 

A. 

where  X(k)  is  the  best  estimate  of  the  states  based  on  past 

information . 

Z(k  -  1)  =  H  X(k  -  1)  +  V(k  -  1)  (19) 

Z  is  the  noisy  observable  matrix,  and  V  is  the  measurement 

noise  vector  defined  by: 

E  [V  ]  =  0  (20) 

E  [Wt>  R  (21) 

R  is  the  covariance  matrix  of  the  measurement  noise  and  P  is 

the  covariance  matrix  of  the  error.    P  is  a  symmetric  matrix, 

(p..   =  p..  ).    The  diagonal  terms  of  P  are  equal  to  the  variance 

of  each  state  (cr^     )  and  the  off-diagonal  terms  are  equal  to 
xi 

correlation  between  the  states, 

(E  [  xi  .  x.  ]-  x.    .  x.). 

P  =  E[  [X-  X].  [X  -  X]1]  (22) 
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The  expected  value  of  the  error  is  defined  as  the  loss  function: 
L  =  I    (X  -  X)t(X  -  X)F(X/Z,X)dX  (23) 

Taking  the  gradiant  of  equation  (23)  with  respect  to  X; 

Vx  L  =    l2XF(X/Z/X)dX  -  2X    J  F(X/Z,X)dX  (24) 

Applying  equations  (16)  and  (17)  and  setting  the  result  equal 

to  zero 

V#L  =  2    E    X/Z,X    -    2  X=  0 

Thus: 

X*  =  Max(x)  =  E[X/Z,X]  (25) 

X*  is  the  best  estimate  of  the  states  given  the  present  noisy 
observable  and  the  past  prediction  of  the  states . 
Filter  Equations 

There  are  two  methods  of  finding  the  recursive  relations  for 
the  filter  program.    Method  one  assumes  the  random  varables  are 
gaussian  and  makes  use  of  Bayes  equation  to  derive  an  ex- 
pression for  equation  (25).    Method  two  assumes  a  linear  fil- 
ter, selects  an  algorithm  for  equation  (25)  and  proves  that  this 
selection  provides  the  best  values  for  X*. 
Method  One 
The  covariance  matrix  of  the  observable  error  may  be  written: 

Pz    =    E[(Z-  ZMZ-zV]  (26) 

Combining  equations  (19)  and  (26), 

Pz   =    E  [  (HX  +  V  -  HX)(HX  +  V  -  HX)*] 
=    H  E  [  (X  -  X)  (X  -  X)1  ]  +  E  [  W*  ] 
Using  equations  (20),  (21)  and  (22), 

Pz    =  HPHt  +  r  (27) 
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In  order  to  find  an  expression  for  X*,  Bayes  formula  is 

applied  to  equation  (25). 

X*  =  F(X/Z/X)  =  F(Z,X/X)  F(X) 

F(Z,X) 

A 
Since  Z(k)  is  independent  of  X(k) ,  equation  (28)  may  be 


(28) 


written: 

F(X/Z,X)  =  F(Z/X)  F(X/X)  F(X)  ,     . 

F(Z)  F(X>  Uyj 

A  *  A 

Since  F(X/X)  =  F(X/X)  F(X)      ,  equation  (29)  becomes, 
F(X) 

F(X/Z,X)  =F(Z/X)  F(X/X)  ,     . 

F(Z)  UU' 

A 

Since  X(k)  is  independent  of  X(k) ,   (30)  becomes, 

F(X/Z,X)  =  F(Z/X)  F(X)  ,     . 

F(Z)  [6i) 

Following  the  form  in  equation  (18),  F(Z)  and  F(Z/X)  are 

defined: 

F(Z)=  1  e-KZ-Z^Pz'^Z-Z)  (32) 

(27>)n/z|Pz|* 

F  (Z/X)  =  F (V)  =  1  e"^vt  R_1  V>  (33) 

(Z$r.)n/2    IRI   * 

Combining  equations  (18),   (31)  and  (33), 

F(X/Z,X)  =A  e"*B 

Where 

A  =         (HPH11  +  Rji 

(2fl)n/^   IR|  8    |P|  * 

and 

B  =  VtR~1V  -  (X  -  X)tP"1(X  -  X)  +  (Z  -  Z^HPH1  +  R)_1(Z  -  Z) 
Setting  the  gradiant  of  B  with  respect  to  X  equal  to  zero. 

VXB  =  ^(X^X)1*  _1  -  2(Z-Z)t(HPHt+R)"1,7x(Z-Z)  =  0 
since 

A 

V    (Z-Z)  =7X(HX  +  V)  =  H,  the  algorithm,  equation  (35) 
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may  be  written. 

X*  =  X  +  PH^HPH1  +  R)"1  (Z  -  Z)  (35) 

Let  the  filter  gain  matrix  G  be  defined. 

G  =  PHMhPH*  +  R)_1  (36) 

Then  the  algorithm,  equation  (35)  may  be  written: 

X*=X+G(Z-Z)  (37) 

Method  Two 

An  alternate  method  of  deriving  the  recursive  relation  for 
the  gain  matrix  assumes  a  linear  filter.    The  trace  of  the  co- 
variance  matrix  is  given  by: 

L  =  E  [  (X  -  X)l{X-  X)]  (38) 

It  will  be  shown  that  if  L  is  minimized  with  respect  to  the  gain, 
then  the  gain  will  be  given  by  equation  (36). 
Substituting  equation  (37)  into  equation  (38) , 

L  =  E  [  (X-X)  (X-X)11  ]     -    E  [  (X-X)  (Z-Z^G*]  - 

E[G(Z-Z)(X-X)t]  +    ECGtZ-ZMZ-zVG1]  (39) 

Since  Z  =  HX  +  V 
Then    Z  =  E  [  HX  +  V  ]  =    HX 

Assuming  that  E  [(X-X)Vt]is  equal  to  zero,  equation  (39)  may 
be  written: 

L  =  E  [  (X-X)  (X-X)11  ]  -  E  [  (X-X)  (X-X)*  ]  HV    +    GE  [ W1  ] 

-  GH  E  [  (X-X)  (X-X)1  ]  +  GHE  [  (X-X)  (X-5?)tHtCt      (40) 
Substituting  equations  (21)  and  (22)  into  equation  (40), 

L  =  P  -  PHkS1  -  GHP  +  GHPHte1  +  GRG1  (41) 

Taking  the  gradiant  of  L  with  respect  to  G  and  solving  for  G. 

G  =  PH^HPH1  +  R)'1 
Filter  Programs 

Jardine,  Titus,  Demetry  and  Strum  have  designed  programs 
to  calculate  the  filter  gains  and  the  covariance  matrix  of  error. 
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Such  a  program,   (Program  Filter)  is  listed  in  Appendix  one. 
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CHAPTER  IV 
DIGITAL  SIMULATION 

If  a  digital  filter  can  be  obtained,  then  a  method  must  be 
developed  to  evaluate  this  system.    This  chapter  will  describe 
the  simulation  of  the  missile  tracking  system,  the  digital  pre- 
dictor, the  control,  and  the  kinematics.    The  results  of  this 
simulation  for  a  particular  missile  will  be  discussed  and  tra- 
jectories from  computer  runs  will  be  presented  in  Chapter  Five. 
Missile  Guidance 

Figure  4-1  is  a  block  diagram  of  the  digital  simulation.    As 
described  in  previous  chapters,  the  missile  system  tracks  the 
target,  obtaining  an  input  for  the  missile  guidance.    This  signal 
is  proportional  to  the  angular  rate  of  change  of  the  line-of-sight 
angle.    In  this  simulation  it  is  assumed  that  white  noise  is 
added  to  this  signal. 

The  details  of  the  missile  guidance  simulation  are  des- 
cribed in  Appendix  Two.    This  portion  of  the  simulation  de- 
scribes the  tracking,  control  and  dynamics  of  the  missile  as 
a  transfer  function  relating  the  line-of-sight  angular  rate  to  the 
missile  flight  path  angular  rate.    This  transfer  function  is  then 
formed  into  an  F  and  D  matrix  as  described  in  Chapter  Three. 
Program  Phidel  is  employed  to  obtain  the  $and  A  matrixes  for 
the  state  difference  equation. 

The  output  of  the  missile  guidance  (y)  is  considered  the 
observable  for  the  predictor.    Gaussian  measurement  noise 
with  standard  deviation  (av)  is  added  to  the  observable.    The 
output  (y)  is  also  an  input  to  the  kinematics.    Figure  4-2 
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Figure  4-1.  Block  Diagram  of  Missile  Simulation. 
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Figure  4-2.  Discrete  Flow  Graph  of  the  Plant. 
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describes  the  digital  simulation  of  the  guidance.    The  following 
fortran  listing  updates  the  discrete  difference  equation  as  dem- 
onstrated in  Figure  4-2  . 
C       THIS  SECTION  UPDATES  THE  STATE  VECTOR  X 

CALL  RNDEV(NUNIF/DEV) 

W=SIGW*DEV 

CALL  PROD  (PHI  ,X,KN,KN,1 , TEMPI) 

DO  803  I  =  1,KN 
803   TEMP2(I,1)  =  W*DEL(I,1) 

CALL  ADD(XS,DINP,KN,1,TELP) 

CALL  PROD(AT,TELP,l,KN,l,TELPl) 

CALL  PROD(DEL  ,TELP1  ,KN ,  1  , 1  ,TELP2) 

CALL  ADD(TEMP1  ,TEMP2  ,KN,  1  ,X) 

CALL  ADD (X,TELP2,KN,1, X) 

Digital  filter-predictor 

The  digital  filter  is  similar  to  the  filter  described  in 
reference    [4  ].    The  gain  matrix  (G)  is  evaluated  each  sample 
instead  of  using  the  steady-state  values  of  the  gain.    The  dis- 
crete flow  graph  of  the  filter  is  shown  in  Figure  4-3. 

One  input  to  the  filter  is  the  noisy  observable.    This  is 
the  sum  of  the  angular  rate  of  the  missile  flight  path  (y)  and 
the  measurement  noise.    This  noise  is  assumed  gaussian  with 
mean  zero  and  variance  ov.    This  measurement  on  the  missile 
could  be  made  by  use  of  rate  gyros . 

The  other  input  to  the  filter  is  the  deterministic  input. 
Here  this  would  be  the  angular  rate  of  the  line-of- sight  (a) 
and  the  best  estimate  of  its  derivatives.    In  the  simulation 
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Figure  4-3.  Discrete  Flow  Graph  of  the  Filter/predictor 
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these  are  obtained  from  the  kinematics.    In  the  missile  this 

signal  would  be  measured  by  the  tracking  system.    Note  that  the 

present  value  of  the  deterministic  input  DI(k)  is  used  as  the 

input  to  the  plant  while  the  past  value  of  the  input  DI(k  -  1)  is 

used  in  the  filter.    The  reasoning  is  that  the  filter  is  taking  the 

present  value  of  the  observable  Z(k)  and  therefore  must  use  the 

past  value  of  the  deterministic  input  DI(k  -  1)  which  caused  the 

observable  Z(k) . 

The  fortran  statements  for  the  filter  are  listed  below. 

C       THIS  SECTION  CALCULATES  XS,  THE  BEST  ESTIMATE 

C       OF  THE  STATE  VECTOR 

CALL  PROD(H,X,KP,KN,l,Y) 

DO  10  1=1, KP 

CALL  RNDEV(NUNIF,DEV) 

10  V(I,l)  =  SIGV(I)*DEV 
CALL  ADD(Y,V,KP,1,Z) 

CALL  PROD(PHI,XS,KN,KN,l, TEMPI) 
DO  11  I  =  1,KN 
DO  11  J  =  1,KN 

11  TEMP2(IJ)  =  -TEMP2(IJ) 

CALL  PROD(TEMP2, TEMPI, KN,KN,1,TEMP3) 
CALL  ADD(TEMP1  ,TEMP3  ,KN  ,  1  ,TEMP3) 
CALL  ADD(XS ,DINP ,KN ,  1  ,TELP) 
CALL  PROD(AT,TELP,l,KN,l  ,TEMP1) 
CALL  PROD(DEL, TEMPI, KN,1,1,TELP) 
CALL  PROD(TEMP2,TELP,KN,KN,l,TELPl) 
CALL  ADD(TELP,TELP1  ,KN,  1  ,TELP1) 
CALL  PROD(G,Z/KN,KP,l/TELP2) 
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CALL  ADD(TEMP3,TELP1  ,KN,1  ,XS) 

CALL  ADD(XS ,  TELP2 ,  KN ,  1 ,  XS) 
The  Kinematics 

The  kinematics  combines  the  flight  path  angular  rate  of  the 
missile,  the  missile  speed  (assumed  constant)  and  the  target's 
velocity  vector  to  determine  their  effect  on  the  line-of-sight 
angular  rate  (a)  and  the  range  rate  (R) . 

Figure  4-4  demonstrates  the  sign  conventions  used  in  this 
simulation.    The  trajectories  which  will  be  discussed  later  in 
Chapter  Five  use  this  same  notation  and  sign  convention. 

Figure  4-5  demonstrates  in  block  diagram  form  the  kine- 
matics.   As  noted  in  the  figure,  by-products  of  the  kinematics 
are  (1)  the  position  of  the  missile  and  target,  (2)  the  range  and 
(3)  the  sign  and  magnitude  of  the  range  rate. 

By  addition  of  a  few  fortran  statements,  adjusting  the 
target's  velocity  vector,  the  target  can  be  maneuvered  as  a 
function  of  the  missile  trajectory  in  range  (evasion)  or  on  a 
predetermined  course  (attack). 

The  fortran  statements  for  the  kinematics  are  listed  below. 
C       THIS  SECTION  DETERMINES  THE  EFFECT  OF  THE 
C       PLANT  OF  THE  KINEMATICS 

GAMDOT  =  X(l,l) 

GAMMA  =  GAMMA+GAMDOT*DT 

XMDOT  =  VM*COSF  (GAMMA) 

YMDOT  =  VM*SLN,  F  (GAMMA) 

YRDOT  =  YTDOT-YMDOT 

XRDOT  =  YTDOT-YMDOT 

RDOT  -  ((YRDOT*SINF  (SIGMA)  )+(XRDOT*COSF  (SIGMA))) 

RDSIG  =  ((YRDOT*COSF(SIGMA))-(XRDOT*SINF  (SIGMA))) 
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Figure  4-4.  Geometry  of  the  Kinematics 
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Figure  4-5.    Block  diagram  of  the  Kinematics. 


RLOS  =  RLOS+RDOT*DT 

DSIG  =  RDSIG/RLOS 

SIGMA  =  SIGMA+DSIG*DT 

XMZ  =  XMZ+XMDOT*DT 

YMZ  =  YMZ+YMDOT*DT 

XTZ  =  XTZ+XTDOT*DT 

YTZ  =  YTZ+YTDOT*DT 
The  Control 

In  this  simulation  a  control  (AT)  was  selected  through  the 
use  of  program  "OPCON"  which  minimizes  the  states  X(k),  ie, 
"bang-bang  control"  .    Other  forms  of  control  are  also  avail- 
able such  as  placing  a  limit  of  the  fuel  or  energy  expended, 
placing  a  limit  on  the  amount  of  acceleration  allowed ,  or  plac- 
ing a  limit  on  the  magnitude  of  the  angle  and  angular  rate  that 
can  be  measured.    These  options  are  recommended  for  future 
study  and  would,  of  course,  be  necessary  in  making  a  more 
complete  study  of  the  missile  systems. 

For  this  simulation  the  weighting  vector  (A*-)  becomes  a 
constant  vector  which  is  imposed  on  the  difference  between 
the  best  prediction  of  the  states  X*(k),  and  the  deterministic 
input,  DINP(k)  to  determine  the  control  to  be  applied. 

Figure  4-6  is  a  discrete  flow  graph  of  the  entire  simulation. 
The  fortran  listing  of  the  entire  program  may  be  found  in  Appendix 
One. 
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Figure  4-6.    Discrete  flow  graph  of  missile  simulation. 
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CHAPTER  V 
SIMULATION  RESULTS 

In  making  a  study  such  as  the  evaluation  of  a  missile  guid- 
ance ,  it  is  quite  difficult  to  decide  which  parameters  to  vary  and 
which  to  hold  constant.    It  is  even  more  difficult  to  decide  which 
curves  hold  the  most  meaning  as  a  measure  of  success.    Since 
the  ultimate  goal  is  to  hit  the  target ,  trajectories  and  confirming 
print-out  were  used  as  a  primary  measure.    Once  a  hit  was  ob- 
tained with  a  set  of  noise  variances,  the  target  trajectory  was 
varied  to  see  the  effect  on  the  missile.    The  noise  was  also 
varied  to  study  the  amount  of  noise  that  could  be  tolerated. 

The  measurement  noise  was  found  to  be  the  most  critical 
value.    A  standard  deviation  of  measurement  noise  of  0.1  radians/ 
sec  was  found  to  be  the  upper  limit.    Above  this  value  the 
missile  had  some  control,  but  the  missile  trajectories  were  far 
from  desirable.    Operation  near  the  limit  of  measurement  noise 
caused  large  ambiguities  of  the  target  information  at  the  be- 
ginning of  the  flight.    The  target  information  improved  as  the 
range  decreased,  however,  and  in  most  cases  the  missile  was 
able  to  capture  the  target. 

Curves  of  the  normal  acceleration  of  the  missile,  the  filter 
gains,  the  range,  and  the  control  versus  time  were  considered 
and  are  shown  in  the  results  where  these  values  appear  to  be 
critical  or  of  interest. 

This  simulation  assumed  that  the  missile  propulsion  has 
already  burned  out  and  the  missile  is  traveling  at  a  constant 
velocity  (Vm) .    This  simulation  considers  only  the  coplanar 
situation.    The  program  could  be  easily  expanded  to  two  channels 
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and  crosscoupling  could  be  considered.    The  trajectories  may  be 
considered  as  yaw  or  pitch  maneuvers.    If  they  are  considered 
pitch  maneuvers ,  the  effects  of  altitude  variations  and  gravity 
must  be  considered.    Otherwise  it  is  assumed  that  the  missile 
has  the  same  dynamic  reaction  to  a  change  in  either  plane.    The 
gravity  consideration  was  not  considered  in  this  simulation. 

The  target  was  made  to  under-go  several  maneuvers:    (1)  a 
turn,  (2)  a  straight  line  course,  and  (3)  a  left/right/left  turn 
(evasive  course).    Figure  5-1  is  a  graphic  representation  of  the 
missile  and  target  initial  conditions  where  A,  B,  C,  D  represents 
the  target  velocity  vector  for  the  four  cases.    The  numerical  values 
of  all  the  parameters  involved  are  listed  in  Table  5-1 .    The  diff- 
erent maneuvers  of  the  target  are  shown  in  Figure  5-2.    A  par- 
ticular computer  run  will  be  designated  (A- 1-2-0. 1-0.1).    The 
"A"  implies  the  initial  conditions  of  the  target  (position  and 
velocity  direction) .    The  first  *'  1 "  implies  that  the  magnitude 
of  the  target's  velocity  is  1000  ft/sec.    The  "2"  implies  that 
the  target  is  flying  a  straight  line  course  (maneuver  two) .    The 
last  two  numbers  are  the  standard  deviation  of  the  excitation 
noise  and  the  measurement  noise  respectively.    Note  that  all 
angles  in  this  program  are  given  in  radians . 

Figure  5-3  to  Figure  5-2'6  are  presented  below  to  display 
the  filtered  missile  capabilities.    Listed  on  each  figure  are  the 
range  at  which  the  program  was  terminated  (R  <  100  ft.)  and  the 
predicted  miss  distance. 
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APPENDIX  I 


THROUGHOUT  THIS  THESIS  COMPUTER  PROGRAMS  WERE  USED»  MADE 
REFERENCE  TO  OR  EXPLAINED.   THIS  APPENDIX  WILL  LIST- THESE  PRO- 
GRAMS AND  GIVE  A  BRIEF  EXPLANATION  OF  THE  DATA  REQUIRED. 

PROGRAM  FILTER. 

THIS  PROGRAM  CALCULATES  THE  OPTIMUM  STEADY-STATE  GAINS 
G,  FOR  THE  DIGITAL  FILTER.   THE  DATA  CARDS  ARE  EXPLAINED  IN 
THE  INITIAL  COMMENT  CARDS  OF  THE  PROGRAM. 

PROGRAM  FILTER 

C      Dl  ORDER  OF  SYSTEM  IN  12  FORMAT 

C      D2  SAMPLING  INTERVAL  IN  8F10.0  FORMAT 

C      D3  F  MATRIX  BY  ROWS  IN  8F10.0  FORMAT 

C      D4  D  MATRIX  BY  COLUMN  IN  8F10.0  FORMAT 

C      D5   VARIANCE  OF  EXCITATION  NOISE  SIGWSQ  IN  8F10.6  FORMAT 

C      D6  R  COVARIANCE  MATRIX  OF  MEASUREMENT  NOISE  IN  8F10.6  FORMAT 

C       PHI  SYSTEM  TRANSITION  MATRIX 

C       DEL  DISTRIBUTION  MATRIX 

C       G   OPTIMUM  GAIN  MATRIX 

C       H   OBSERVABLE  MATRIX 

C       P  BEST  ESTIMATE  OF  ERROR  COVARIANCE  MATRIX 

C       Q   EXCITATION  NOISE  COVARIANCE  MATRIX 

DIMENSION  P(12»12)»Q(12»12),H(12»12)»R(12»12)»G(12»12)»PHIT(12»12 

It  PHI ( 12tl2) »DEL(12 ),DELDELT( 12»12) »DELS( 12.12) »DELST ( 12 . 12 ) » 

2X(  12)  ,XHAT*(  12)  .YHAT(  12  )  »PNEW(  12*12) 

READ  33.  N 

33  FORMAT  (12) 

READ  2,DT 

2  FORMAT(8F10.0) 

PRINT  1003 

1003  FORMATt 1H1) 

CALL  PHIDEL(PHI ♦DEL.N.DT) 

DO  1001  LL=1.A 

READ  20.SIGWSO 

20  FORMAT(F10.6) 

DO  1001  LK=1,4 

C      INIALIZE  MATRICES  FOR  EACH  RUN 

DO  10  I=1,N 

DO  10  J=1»N 

PNEW( I »J) =0.0 

P(  I  »J)=0.0 

G(  I  »J)=0.0 

0(  I  ^5=0.0 

DELDELT( I, J)=0.0 

0FL6.(  t  Mjso*0 

10  DELST( I »J)=0.0 

READ  2001,R(1»1 ) 
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PROGRAM  FILTER  CONTINUED 

2001  FORMAT(8F10.6) 
00  30  I»1»N 

30  delsu»1)=del<  t> 

CALL  TRANS(DELStNtNtDELST) 
CALL  PROD(DFLS.DFLST»N«N.N.DELDFLT> 
00  40  I=1»N 
DO  40  J  =  1#N 
40  Q< I.J)=SIGWSQ*DELDELT< I»J) 
SIGW=SQRTF(SIGWSQ) 
PRINT  1004.SIGW 

1004  FORMATf/  10Xt5HSIGW=  /»E20.8) 
P(ltl)=.02 

P(lt2)=0.0 
P(2»l)=0.0 
P(2»2)»1.0 
PRINT  2002.R(1.1W<  ( Q(  I  •  J )  »  J*l  .N)  1 1 *1  »N)  . ( (P(  I  •  J)  . J»l  tN)  •  I*  1  »N) 

2002  FORMAT(///20X»8HR(1,1)*  »F10.6/20X t8H0( I ♦ J  J =  •4F10«6/20Xt8HP( ! t J)« 
1  .4F10.6) 

H(ltl)«l« 

H(  l#2)s»0.0 

PRINT    700 
700    FORMAT (//5X.3HGll,7X»3HG12.7X,3HG21.7X»3HG22#7X»3HPll#7Xt3HP12t7Xt 
1    3HP2U7X.3HP22/) 

DO    1000    KK=lt40 

CALL    GP(H.PHI»P»Q,Rf2»l.G»PNEW) 

PRINT  800 ( (G( I.J)»J»l»N)tI»ltN) •  < (PNEW( I • J) • J«l .N) •I«1*N) 
800  FORMAT  (10F10.5) 

DO  11  I=1»N 

DO  11  J=1»N 
11  P(I»J)*PNEW(ItJ) 

1000  CONTINUE* 
PRINT  1005 

1005  FORMAT(lHl) 

1001  CONTINUE 
END 

• 

SUBROUTINE    GP < HtPHI ♦ P,Q»R tKN.KP tGt PNEW) 

DIMENSION    H(12»12).PHI(12tl2)tP(12.12)tQ(12»12)»R(12.12)»G(12tl2l» 
1PNEW(12»12) •  HTU2.12WTVH12.12)  #TV2<12»12> 
CALL    TRANS(H,KP»KN»HT) 
CALL    PROD(P.HT,KN»KN,KPtTVl) 
CALL    PR0D(H,TV1»KP»KN»KP»TV2) 
CALL    ADD(TV2tR.KP»KP#TVl) 

CALL    REC I P(KP», 0000000000001 »TVl»TV2tKER) 
IF<KER-2)     101»110tl01 

110  PRINT  111 

111  F0RMAT(5HKER*2) 

101    CALL  PR0D(HT»TV2»<N»KP.KP»TV1) 

CALL  PROD(P»TVl»KNtKN»KPtG) 

(CALL  P03d<NtPa|fcPaftNiKN»TVH 

CALL  PR0D(G»TVl.KN.KP»KNtTV2) 

73 


PROGRAM  FILTER  CONTINUED 

DO  102  1  =  1, KN 
DO  102  J=1,KN 
102  TV2( I ,J)=-TV2( I ,J) 

CALL  ADD(P,TV2,KN,KN,TV1) 

CALL  PROD(PHI , TVl » KN , KN , KN , TV2 ) 

CALL  TRANS(PHI ,KN,KN,TV1 ) 

CALL  PR0D(TV2,TV1,KN,KN,KN,PNEW) 

CALL  ADD(PNEW,Q,KN,KN,PNEW) 

END 


153 


SUBROUTINE  TRANS ( A ,N »M ,C ) 

DIMENSION  A(12, 12) ,C( 12,12) 

DO  153  I  =  1,N 

DO  153  J=1,M 

C(  J»I )  =  A( I ,J) 

END 
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2 

10 


11 


12 


13 


14 


15 
20 
30 
31 


32 

33 

35 
36 


SUBROUTINE  REC I P ( N , EP , A ,X ♦ KER ) 

DIMENSION  A(12,12) ,X(12,12) 

DO  1  I=1,N 

DO  1  J=1,N 

X(  I  ,J)=0.0 

DO  2  K=1,N 

X(K,K)=1.0 

DO  34  L=ltN 

KP  =  0 

Z  =  0.0 

DO  12  K  =  L,N 

IF(Z-ABSF(A.(K.L)  )  )  11  » 12  #12 

Z=ABSF(A(K,L) ) 

KP  =  K 

CONTINUE 

IF(L-KP) 13,20,20 

DO  14  J=L,N 

Z=A(L,J) 

A(L,J)=A(KP,J) 

A(KP»J)=Z 

DO  15  J=1,N 

Z=X(L,J) 

X(L,J)=X(KP,J) 

X(KP,J)=Z 

IF(ABSF(A(L,L) ) -EP ) 50 , 50 ♦ 30 

IF(L-N)31»34,34 

LP1=L+1 

DO  36  K=LP1,N 

IF(A(K»L)  )32, 36,32 

RATIO=A(K,L)/A(L»L) 

DO  33  J=LP1,N 

A(K,J)=A(K,J)-RATIO*A(L»J) 

D6  3£  JsiiN 

X(K,J)=X(K, J)-RATIO*X(L,J) 

CONTINUE 
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PROGRAM  FILTER  CONTINUED   " 

34  CONTINUE 

40  DO  43  I=lfN 
II=N+1-I 

DO  43  J  =  1.N 

S  =  0.0 

IFC II-N)41»43»43 

41  IIP1=II+1 

DO  42  K=IIP1»N 

42  S=S+A( II»K)*X(K.J) 

43  X( I I.J)=(X< I  I » J)-S) /A( I  I » 1 1 ) 
KER  =  1 

RETURN 
50  KER=2 
END 


SUBROUTINE  PHIDELt PHI ♦DELtN»DT ) 
VALID  ONLY  FOR  A  CONSTANT  F  MATRIX 

DIMENSION  F( 12*12) tPHI ( 12 • 12 ) ,TERM( 12*12) •WORM! 12*12) 
It  DEL (12) tDELM( 12.12) tTELM( 12 • 12 ) .DELP( 12*12 )*0( 12) 
READ1.( <F<IRtIC)tI01*N)*IR*l*N) 
1  FORMAT  ((8F10.0)) 
READ  1  (D(I )»I*1.N) 
1003  PRINT  399»DT»((FUR*IC)*IC*1*N)»IR-1*N) 

399  FORMAT(///3HDT=r.lF5.3///»7HF(I,J)«/.( (8F8*2) ) ) 
PRINT  3991  (D(I ),I=1»N) 

3991  FORMAT(////5HD( I)*/(8F8.2) ) 
NFINAL*1 
TM=0.0 

DO  400  IR=1,N 
DO  400  IC=1.N 
TERM(IR.IC)=0.0 
WORM(IR»IC)=0.0 
TERM(IR»IR)=1.0 
TELM(IR.IC)=TERM( IR.IC)*DT 
DELP(IR.IC)=TELM(IR.IC) 
DELM( IR.IC)=0.0 
DEL( IR)=0.0 

400  PHI  UR.IC)=TERM(IR»IC) 
4  TM=1.0+TM 

DO  500  IR=1.N 
DO  500  IC=1.N 
DO  500    JN*1»N 

DELM(IR.IC)=DELM(IR»IC)-TELM( I R» JN )*F(JN» IC ) *DT/< TM+1.0 > 
500  WORM(IRtIC)=TERM( IR. JN )*F( JN. IC )*DT/TM+WORM( IR*IC) 
DO  401  IR=1.N 
DO  401  IC=1»N 
TERMUR.IC)=WORM(  IR.IC) 
TELM( IR.IC)=DELM( IR*IC) 
DELP( IR.IC)=DELPUR.IO  +  TELM<!R.IC) 

DELM(IR.IC)=0.0 

401  WORM( IR.IC)=0.0 
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ABC=0.0 

DO  2I=1»N 

DO  2J=1»N 

AA=TERM( I ,J) 

AB=ABSF(AA) 

IF(ABC-AB)3,3*2 
3  ABC=AB 
2  CONTINUE 

IF(0.O00000  00  5-ABC)5,5.6 

5  GO  TO  4 

6  PRINT  502,(  (PHI  (  IR,IC)  »I01»N)  ,IR=1,N) 
50  2  FORMAT <////9X,8HPH  I (  I  » J  ) /// ( 6E20.8 )  ) 

DO  600  I=1»N 

DO  600  K-ltN 

DO  600  J=1»N 
600  DFL( I )=DEL( I )+PHI( I , J ) *DELP ( J , K ) *D ( K ) 

PRINT  503  (DEL( I ) ,I=1.N) 
50  3  FORMAT (////9X,6HDEL(  I ) // ( 6E20.8 ) // ) 

END 


SUBROUTINE  ADD  (A»B»N»M»C) 
DIMENSION  A( 12.12) »B ( 12 .12 ) *C( 12*121 
DO  152  1  =  1, N 
DO  152  J=1.M 
152  C( I »J)  =  A( I ,J)  +  B( I »J) 
END 


SUBROUTINE  PROD  ( A.B.N»M»L »C ) 
DIMENSION  A(12»12)*B(12.12)»C(12»12) 
DO  151  I=1»N 
DO  151  J=1,L 
C(I»J)  =0. 
DO  151  K  =  1»M 
151  C(I»J)  =  C(I.J)  +  A(  I  ,K)*B(K,J) 
END   . 
END 


PROGRAM  OPTCON 
DATA  CARDS 

Dl  CASE  THE  COMPUTER  IS  TO  EXAMINE. 

CASE  1Q=I.R=0.Q*=0 

CASE  2Q=0,  R=1*Q*=0 

CASE  3  Q  =  I.  R  =  1.  Q*  EMPLOYED 
D2  N  =  ORDER  OF  SYSTEM*  NSTAGE  =  NO.  OF  ITERATIONS 
D3  O-WEIGHTING  MATRIX 
D*  F  MATRIX 
D5  D  VECTOR 
D6  DT-SAMPLE  PERIOD 

7$ 


PROGRAM  OPTCON 

DIMENSION  PHI (12 ,12) .PS H 12 ,12) *P( 12*12 ) *DEL ( 12 )• ATI  12) • 
1GM( 12.12) »Q( 12*121 »X(900).ITITLE( 12)*FM< 16 ) *EM( 12 ) *HM( 12* 12 ) • 

1  Y1(900)»Y2(900).Y3(900) 

C  THIS  PROGRAM  UTILIZES  A  COST  FUNCTION.  J(N )  *MINIMUM<  SUM  X< N> T*Q*X( N>* 

C  SUM  R*U(N-1)**2).   AN  UNLIMITED  NUMBER  OF  ITERATIONS  MAY  BE  MADE  AT 

C  A  COMPUTATION  RATE  OF  2000  PER  MINUTE  AFTER  THE  PROGRAM  HAS  BEEN 

C  COMPILED.   THE  OUTPUT  OF  THIS  PROGRAM  IS  THE  FEEDBACK  GAIN  MATRIX* 

C  A  TRANSPOSE.   THE  FOLLOWING  RECURSIVE  EQUATIONS  WERE  DERIVED  USING 

C  DYNAMIC  PROGRAMMING. 

C  AT(K)=-(DELT*P(K-1)*PHI ) / <DELT*P< K-l ) *DEL+R )  (1) 

C  PSI (K)=PHI+DEL*AT(K)  (2) 

C  P(K)=PSIT(K)*P(K-1)*PSI (K)+Q+R*A(K)*AT{K)  13) 

C      P(0)=0.AT(0)=0.  PSI(0)=0 

C  EQUATIONS  1.  2.  AND  3  CONSTITUTE  THIS  PROGRAM 

C      CALL  IN  DATA  AND  INITIALIZE 
DO  1111  1=1.3 
READ  30.KASE 
30  FORMAT ( II) 

READ  1  .N.NSTAGE 

1  FORMAT  (8110) 
READ  2.R.DT 

2  FORMATt (4F20.0) ) 

READ  2  ( (Q( I.J) .J*1*N) *I>1*N) 
C      PRINT  THE  DATA  READ  IN* 
PRINT  100 
100  FORMAT(lHl) 

PRINT  32.KASE 
32  FORMAT (9X.5HKASE*  .13) 
PRINT  3. N.NSTAGE 

3  FORMAT(//9X.2HN=.I3*20X*7HNSTAGE»*I3) 
PRINT  4.R.DT 

4  FORMAT(/9X.2HR=.F15.11.2X.3HDT=»F15*ll) 
PRINT  9.( (Q( I.J).J=1.N).I*1.N) 

9  FORMAT </9X.7HQ( I » J )«/ ( 2F15.11 ) ) 
CALL  PHIDEL(PHI.DEL.N.DT) 
DO  5  I=.1.N 
DO  5  J=1.N 
GM( I.J)=0.0 
EM( I )=0.0 
FM( I )=0.0 

5  P( I.J)=Q( I.J) 

C      INITIALIZE  P(0)  FOR  CASE  TWO 
IFUASE-2)  35.36*35 
36  P(l.l)=1.0 
P(2.2)*1.0 
P(3.3)=1.0 
35  CONTINUE 
PRINT  19 
19  FORMAT ( //2X,6HNSTAGE.4X»7HAT(1.1)*3X*7H AT (1*2)*4X*6HP( 1*1)* 
1*X*6HPU*2)  »AX*6HP(2*1  )  .4X  *6HP  (  2  *2  )  ) 
C      CALCULATE   AT(J) 

DO   22  KK=1.NSTAGE 

?l 


PROGRAM  OPTCON  CONTINUED 

DEN=0.0 
D06  I  =  1»N 
D06  J=1,N 

6  EM( I )=EM( I )+DEL(J)*P( J,I  ) 
D08  1  =  1, N 

D07  J=1,N 

7  FM(  I  )=FM(  I )+EM( J)*PHI ( J, I  ) 

8  DEN  =  DEN+EM(  I )*DEL(  I  ) 
DEN=-DEN-R 
DO10I=l,N 

AT( I )=FM( I )/DEN 

FM( I )=0.0 
10  EM( I )=0.0 
C      CALCULATE  PSKKtltJ) 

D013I=1.N 

D013J=1»N 
13  PS  I ( I »J)=PHI (  I ,J)+DEL(  I  )*AT( J) 
C      CALCULATE  P(K»I»J) 

DO  15  I=1,N 

DO  15  J=1»N 

DO  15  L  =  1,N 

15  GM( I ,J)=GM( I >J)+PSI (L»  I  )*P(L»J) 
DO  171=1, N 

DO  17J=1,N 
DO  16  L=1,N 

16  HM(  I  ,J)=HM(  I  ,J)+GM(  I  ,L)*PSI  (L,J) 

C      CASE  1   TERMINAL  CONTROL,  OMIT  Q(I,J)  IN  EQUATION  FOR  P(I,J) 
C      CASE2   MINIMUM  FUEL        OMIT  Q(I»J)  IN  EQUATION  FOR  PII.JJ 
IF(KASE-2)  31,31,33 
31  P(  I  ,J)=HM(  I  ,J)+R*AT<  I )*AT( J) 
GO  TO  17 

C      CASE  THREE  MINIMIZATION  OF  TIME  AND  FUEL 
3  3  P(  I  ,J)=HM( I ,J)+Q(  I  ,J)+R*AT( I )*AT(J) 

17  HM(  I  ,J)=0.0 
D018  I-=1,N 
D018  J=1,N 

18  GM(  I  ,J)=0.0 

PRINT  20,KK,AT( 1) ,AT(2) ,P( 1,1) ,P( 1,2) ,P(2,1) ,P(2,2) 
20  FORMAT ( I5,(6F10.6) ) 
22  CONTINUE 
1111  CONTINUE 
END 

SUBROUTINE  PH I  DEL ( PH I , DEL ,N , DT ) 
C      VALID  ONLY  FOR  A  CONSTANT  F  MATRIX 

DIMENSION  F( 12, 12), PHI (12, 12) ,TERM( 12,12) ,WORM(12,12) 
1,  DEL(  12) ,DELM( 12,12 ) ,T ELM (12, 12) »DELP( 12,12) ,D( 12) 
1  FORMAT  ((8F10.0)) 

READ1, ( (F(IR,IC),IC=1,N),IR=1,N) 
R^AD  1  (D( I) ♦ I=1»N) 
1003  PRINT  399,DT,( (F( IR, IC) ,IC=1,N) ,IR=1,N) 
399  F0RMAT(///3HDT=,1F5.3/   , 7HF ( I , J ) =/ , ( ( 2F8.2 ) ) ) 
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PRINT  3991  (D( I ) .I=1.N) 
3991  FORMAT(/    5HD( I ) */ < 2F8.2 >  ) 
NFINAL=1 
TM=0.0 

DO  400  IR=1.N 
DO  400  IC=1.N 
TERM(IR.IC)=0.0 
WORM(IR»IC)=0.0 
TERM( IRtIR)=1.0 
TELM( IR.IC)=TERM( IR. IC)*DT 
DFLP( IR.IC)=TELM( IR.IC) 
DELM( IR.IC)=0.0 
DEL( IR)=0.0 

400  PHI  <  IR.IC)=TERMUR.IC) 

4  TM=1.0+TM 

DO  500  IR=1»N 
DO  500  IC=1»N 
DO  500    JN=1»N 

DELM( IR.IC)=DELM( IR. IC)-TELM( I R. JN) *F< JN. IC >*DT/< TM+1.0 ) 
500  WORM( IR,IC)=TERM{ I R. JN )*F { JN. IC )*DT/TM+WORM< IR.IC) 
DO  401  IR=1,N 
DO  401  IC=1.N 
TERM( IR.IC) =WORM( IR.IC) 
TELMt IR»IC)=DELM(IR» IC) 
DELP(IR.IC)=DELP(IR.IC)+TELM< IR.IC) 
PHI( IR»IC)=PHI( IR.IC)+TERM( IR.IC) 
DELMUR,  10=0.0 

401  WORM( IR»IC)=0.0 
ABC=0.0 

DO  2I=1»N 

DO  2J=1.N* 

AA=TERM( I.J) 

AB=ABSF(AA) 

IF(ABC-AB)3,3»2 
3  ABC=A8 
2  CONTINUE 

IF(0.000000005-ABC)5.5.6 

5  GO  TO  4 

6  PRINT  ll.( (PHI( I »J)» J=1.N) .I=1.N) 

11  FORMAT (/9X.9HPH I <  I»J)=/(2F15.11)) 
DO  600  1=1. N 

DO  600  K=1.N 
DO  600  J=1.N 
600  DEL(  I  )=DFL(  I  )+PHI(  I .  J  )  *DELP  (  J,K  )  *D(  K  ) 
PRINT  12.(DELU ).I=1.N) 

12  FORMAT (9X.7HDEL( I )=/(l F15.ll >) 
END 

END 
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PROGRAM  MISSILE 

C  THIS  PROGRAM  SIMULATES  THE  MISSILE  CONTROL ,DYNAM I CS > AND  KINEMATI' 

THE   KALMAN  FILTER  AND  OPTIMAL  CONTROL  IS  APPlIED  TO  THE  CONTROL 

C  OF  THE  PLANT. 

C  THIS  PROGRAM  USES  VARIABLE  FILTER  GAINS 
C 

C  THIS  PROGRAM  CLOSES  THE  i_OGP  ON  THE  OPTIMUM  FILTER-CONTROLLER 

C  PROBLEMS.  IT  ASSUMES  THAT  STABLE  CONTROLLER  GAIN  MATRIX  HAS 

C  BEEN  COMPUTED  ON  THE  BASIS  OF  DESIRED  RESPONSE  AND  CALCULATES 

C  THE  FILTER  GAIN  MATRIX  EACH  SAMPLE  ON  THE  BASIS 

C  OF  THE  STATISTICAL  PROPERTIES  OF  THE  ANTICIPATED  RANDOM  DISTURBAI 
C 

C  THE  PROGRAM  SOLVES  THE  FOLLOWING  EQUATIONS 
C  Y(K)=H*X(K) 

C  Z(K)=Y(K)+V(K) 

C  XS(K)=(  I-GH)PHI*XS(.<-1  )  +  {  I-GH)-»-DEL*AT(XS(K-l  ) -DINP(  K-l )  )+G*Z  (  K)' 

C  X(K+l)=PHI*X(K)+DEL*ta ( K)+DEL*AT*(XS( K)-DINP« K) ) 

C  WHERE  V  IS  MEASUREMENT  NOISE, W  IS  THE  RANDOM  DISTURBANCE » AND 

C  DINP  IS  THE  DETERMINISTIC  INPUT 
C 

DATA  CARDS. 

Dl  N  =  ORDER  OF  SYSTEM,  DT=  SAMPLE  PERIOD 
D2  GAMMA  =  FLIGHT  PATH  ANGLE 

SIGMA  =  LOS  ANGLE 

GAMDOT  =  FLIGHT  ■      TH  ANGULAR  .'RATE 

DSIG  =  LOS  ANGULAR  RATE 

D3  RLOS  =  RANGE 

RZ  =  REQUIRED  TARGET  RANGE  FOR  HIT 

XTZ  =  INITIAL  X-COORDINATE  OF  TARGET 

XTDOT  =  INITIAL  VELOCITY  IN  X  DIRECTION  FOR  THE  TARGET 

YTZ  =  INITIAL  Y-COORDINATE  OF  TARGET 

YTDOT  =  INITIAL  VELOCITY  IN  Y  DIRECTION  FOR  THE  TARGET 
DA  INITIAL  CONDITIONS  OF  THE  STATE  VECTOR 
D5  XMZ  INITIAL  X-COORDINATE  OF  MISSILE 

YMZ  INITIAL  Y-COORDINATE  OF  MISSILE 

VM  MISSILE  SPEED 
D6  P-BEST  ESTIMATE  OF  ERROR 
D7  F  MATRIX 
D8  D  VECTOR 

D9  GRAPH  TITL.E  LINE  ONE 
D10  AT  FEEDBACK  GAIN  MATRIX  FOR  CONTROLLER 

.  1  GRAPH  TITLE  LINE  TWO 
D12  SIGW  =  STANDARD  DEVIATION  OF  EXCITATION  NOISE 

SIGV  =  STANDARD  DEVIATION  OF  MEASUREMENT  NOISE 

D  I  MENS  I  ON  X  (  12  » 12  )  ,XS  ( 12  » 12  )  » S  IGV  (  12  }  » Y  (  12..  12  )  *Z  (  12  »  12  )  » PH  I  (  12  ,  12 
1»DEL(12»12)>H(12>12)»AT(12»12),G(12,12)  ,TEMP1( 12,12)  »TEMP2(  12,12) 
2TEMP3(  12»12)  ,TEMPA( 12,12),TELP(12,12)»TELP1{12,12),TELP2(12,12)» 
3V(  12,12!  ,DINP(  12,12)  ,  IT ( 12) ,TME( 900)  ,YT (900)  ,  XT (900) 
A,XM(900) ,YM(900 ),P(i2,12),R(l2,12) 
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PROGRAM  MISSILE  CC  T   JED 

5»DELTR( 12,12) »DELS( 12) »QC 12»12) 

6»SDOT(900} »DIFF(900) 
C      INPUT  SOME  CONSTANTS*  AS  NOTE D 

READ  2,N,DT 

READ  1 »GAMMA»SIGMA»GAMDOT»DSIG 

READ  l»RLOS»RZ,XTZ»XTDQT»YTZ»YTDOT 

READ  1»(X(I»1)»I=1»N) 

READ  1»XMZ»YMZ,VM 
3  FORMAT (6A8) 
2  FORMAT  (  15  -,  1F10.0) 
1  FORMAT ( (8F10.0) 5 

READ  1»<{P(I»J)»U=1, N 5  ,  I  =  1  ,  N  ] 

CALL  PHIDEUPHI  »DELS»N*DT) 

READ  3, (ITU), 1=1, 61 

READ  1  * (AT  £ 1» I )  ,1  =  1  ,H) 

READ  3»( IT( I ) >  1=7*12) 

READ  1»SIGW»SIGV( 1 ) 

PRINT  144 
144  FORMAT(4X>8HXMISSILE»6X»8HXTARGET  i  8X»  8HYMISSILE»8X  »8HYTARGET  ) 

DO  6  00  0  I=1,N 
6000  DFL(  I  »  1  )=DFLS( I  ) 

CALL  TRANS(DEL,N,N,DELTR] 

CALL  PROD(DEL»DELTR,N  »  N*N>Q) 

DO  2  00  J=] ,N 

DO  2  00  1=1, N 
20  0  Q{ I ,J)=Q( I , J)*SIGW*SIGW 

R( 1,1)=SIGV£1)*SIGV{  I ) 

P(1»1)=R£1,1) 
C      :<N  IS  SYSTEM  ORDER,  KP  IS  THE    NUMBER  OF  OBSERVABLES. 

KN  =N 

KP  =  1 

H(l»l)=1.0 
C      SETTING  DINP=0.0  AT  "r.lS    POINT  INSURES  THAT  NO  DETERMINISTIC 
C      INPUTS  EXIST  PRIOR  TO  TIME  =  ZERO 

DO  131  1=1, KN 
131  XS( I »1 )=0.0 

NUNIF=1220703'125 
C      FILTER  INITIAL  CONDITIONS 

C-n-L  PROD(H,X,KPsKN,KP,Y) 

DO  12  1=1, KP 

CAuL  RNDEV  (  NUN  I  F  •,  D  EV  ) 
12  V<  I  >1)=SIGV( I )*DEV 

CALL  ADD(Y»V,KP,i,Z) 

XS  C 1 » 1 5 -Z (1,1) 

KK=1 

T=0.0 

KK1=0 
C      THIS  SECTION  CALCULATES  G,  THE  GAIN  MATRIX 

600  9  CALL  F I LTER  (  KN ».<P»PH I  »Q»R»H,P,G) 

C 

C      THIS  SECTION  CALCULATES  XS»  THE  BEST  ESTIMATE 
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10 
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VECTOR 
,KP,KN, 1,Y) 

NIF,DEV) 
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KP, 1,Z) 
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,KN,KP,KN»TEMP2) 
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1  .. 
LP 
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PI 
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■. 
KNj 
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N»KN,1 » TEMP 3) 
,  1,TEMP3) 
ELP) 

1»TEMP1 j 
1»1"»TELP) 
,KN,1,TELP1 ) 
1 »TELPl ) 
Z_P2  ) 

»  1  9  X  S  ) 

XS) 


DINP( 1 


DSIG 


DIMP( 2  » I) =-DDSIG 

SDOT(KK)=DSIG 

TME(KK)=T 

THIS  SECTION  UPDATES 


STAIE  Vi 


OR  X 


CALL  RNDEV<NUNIF,DEV) 

W=SIGW*DEV 

CALL  PROD ( PH I » X  » KN , KN , 1 , T EMP  \ 


80 


TEMP2( I .1 )=W*DEL( I  si ) 
CALL  ADD(XS,DINP,KN»1»TELP) 
CALL  PROD ( AT  ,TELP, 1 » KN , 1 »TELP1 ) 
CALL  PROD (DEL  » TELP1 , KN » 1 » 1 ,TELP2) 
CALL  ADD ( TEMP 1,TEMP2 ,KN,1,X) 
CALL  ADD(X,TELP2,KN,1»X) 


THIS  SECTION  DETERMINES 

ON  THE  KINEMATICS 
DKSIG=DSIG 


THE  EFFECT 


He. 


.ANT 


DIFF(KK)=X( 1 »1)-XS (1,1 ) 

GAMDOT=X( 1,1 ) 

GDOT(KK) =GAMDOT 

GAMMA    =GAM  IA+G  ■     Z  DT*DT 

XMDOT=VM#COSF(G- 
13  3  YMDOT=VM*S I NF ( GAMMA ) 
134  YRDOT  =YTDOT-YMDOT 

XRDOT  =XTDOT-XMDOT 

RDOT=( (YRDOT*S INF (SIGMA) ) + ( XRDOT*COSF(S I GMA ) ) ) 

RDSIG=( (YRDOT*COSF< SIGMA) )-£ XRDOT*SINF( SIGMA ) ) ) 


PROGRAM  MISSILE  CONTINUED 

RLOS  =RLOS+RDOT*DT 
DSIG=RDSI    JJ$ 
SIGMA  =SIGMA  -DSIG*DT 
DDSIG=(DSIG-DKSIG)/DT 
a  4Z=XMZ+XMDOT*DT 


140 
141 


190 

145 

146 

6012 

6011 

6013 

13  5 

136 


XM(KK)  -XMZ 
YM(KK) =YM2 
XTZ  =XTZ+XTDOT*DT 
YTZ  =YTZ+YTDOT*DT 
YT(KK)  =YTZ 

xt  ::<;<}  =xtz 
rl(kk)=rlos 

PRINT  140 »KK 

PRINT  ,  1 41,  XM  {,<;<)  ,XT 

FORMAT(2X,I5) 

FORMATt  2X,4E16.7) 

KK=KK+1 

T  =  T  +  DT 

PRINT  190,  DSIG*SIG 

FORMAT(6H  DSIG=*1E1 

IF (ROOT) 146 , 1  h5 .145 

PRINT  135»RDOT 

GOTO  6013 

IF(KK-900  )6012»6012,6013 


KK)  >YMCKK)  ,YTCKK) 


2X»7H  SIGMA' 


:i6.7i 


GOTO  6009 
PRINT  1 3 6  t 


RL< 


FOR  M  A  i  (  6  H  RD0TS!>lE16e7) 
F0RMAK6H  R.OS^  ,  IE  16.  7  J 
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151 


TMISS=(  CRLOS*RDSIG)/RDOT)*COSF( A  J 

PRINT  137,TMISS 

FOR  AT(7H  TMISS-*iE16»7) 

1/  U  _.  '   ' 

N  N  -  l\  .N  —  1 

MC  -1 

LA=  LH  ,-,  \ 

CALL  DRAW{  KK»XM»YM»MC  ?0»LA»  IT,0 

MC-3 


»J*vv»U»0 


C A  _  _  0 RAW  v  KK  v  X  i  » Yl  -j  MC  »  0  »  LA  » I T  *  0 

E  M  D 

SUBROUTINE  PROD  (A*3»N»M*L»C) 

D [ME  45 1  ON  A( 12*12)  » 3 { 12 » 12) »C( 12*12 } 


0»0,0»Q»0»7*S»0,LAST  J 
>  »G  *7»  S  »  G  »  LAST  ) 


DO    151     I=1,N 

DC    151    J=1»L 

DO    151    K=1*M 

C(I»U)=:C(I*«J/"«"  A 

;  i ,  k  ) 

*B(K»J 

END 

S U  b K  OU  ,  I N  E    A  -/  -' 

i  A  v  b  j 

N  » ivU  C  ) 

DIMtNS . ON 

DO  152  „- 


2*12) »C£ 12,12) 
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DO  152  1=1, N 
152  C( I »J)=A( I . J)+B( I ♦ J) 
END 


SUBROUTINE  F I LTER ( N, KP  ♦  PHI ,0 , R ,H ,P ,G ) 
C       PHI  SYSTEM  TRANSITION  MATRIX 
C       DEL  DISTRIBUTION  MATRIX 
C       G   OPTIMUM  GAIN  MATRIX 
C       H   OBSERVABLE  MATRIX 

C       P  BEST  ESTIMATE  OF  ERROR  COVARIANCE  MATRIX 
C       Q   EXCITATION  NOISE  COVARIANCE  MATRIX 

DIMENSION  P( 12*12)  *Q(  12.12  )-*H(  12.12)  »R(  12  . 12)  »G(  12  » 12)  .PHI  T(  12*] 
1  .PHI  ( 12.12)  ,DEL( 12) , DELDELT ( 1 2 , 12 ) ,DELS( 12.12) .DELST ( 12  ♦  12  )  . 
2PNEW( 12.12) 

CALL  GP(H.PHI»P.O.R.M»KP,G.PNEW) 

DO  11  I=1»N 

DO  11  J=1,N 
11  P( I .J)=PNEW( I.J) 

END 

SUBROUTINE  GP ( H , PH I , P , 0 ,R »KN . KP .G . PNEW ) 
DIMENSION  H(12,12).PHI(12»12),P(12.12).0(12.12).R(12.12).G(12.12) 
1PNEW<12»12) »HT( 12.12) »TV1( 12.12 ) »TV2< 12, 12) 

CALL  TRANS(H.KP.KN.HT) 

CALL  PR0D(P,HT»KN»KN.KP»TV1) 

CALL  PR0D(H,TV1,KP,KN,KP,TV2) 

CALL  ADD(TV2»R»KP.KP»TV1) 

CALL  RECIP( KP.. 0000000000001 .TV1.TV2.KER) 

IF(KER-2)  101,110,101 

110  PRINT  111 

111  FORMAT(5HKER=2) 

101  CALL  PR0D(HT.TV2.KN,KP.KP.TV1) 
CALL  PROD(P.TVl.KN.KN.KP.G) 
CALL  PR0D(H»P,KP,KN»KN.TV1 ) 
CALL  PR0D(G,TV1»KN»KP»KN.TV2) 
DO  102  I=1.KN 
DO  102  J=1,KN 

102  TV2( I ,J)=-TV2( I ,J) 
CALL  ADD(P»TV2»KN»KN.TV1) 
CALL  PROD(PHI . T VI , KN , KN , KN , T V2 ) 
CALL  TRANS(PHI .KN.KN.TV1) 
CALL  PROD (TV2.TV1.KN.KN.KN, PNEW) 
CALL  ADD(PNEW.Q.KN.KN.PNEW) 
END 

SUBROUTINE  TRANS ( A ,N »M ,C ) 
DIMENSION  A( 12,12)  ,C(  12.12) 
DO  153  I  =  l.N 
DO  153  J=1.M 
193  £(  Jtli  =  hi  ;  »~m 
END 
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SUBROUTINE  TRANS ( A ,N .M ,C ) 
DIMENSION  A(  12.12)  .C( 12.12  ) 
DO  153  I  =  l.N 
DO  153  J=1.M 
153  C( Jtl )  =  A( I .J) 
END 

SUBROUTINE  REC I P ( N ,EP , A.X.KER ) 
DIMENSION  AU2.12)  .XU2.12) 
DO  1  1=1 .N 
DO  1  J=1»N 

1  X(  I  ,J)=0.0 
DO  2  K  =  1.N 

2  X(K«K)>i«0 

10  DO  34  L=1.N 
KP  =  0 

Z  =  0.0 

DO  12  K  =  L.N 

IF(Z-ABSF(A<K»U )) 11.12.12 

11  Z=ABSF( A(K.L) ) 
KP  =  K 

12  CONTINUE 
IF(L-KP)13.20,20 

13  DO  14  J=L»N 
Z=A(L.J) 
A(L.J)=A(KP,J) 

14  A(KP,J)=Z 
DO  15  J=1.N 
Z=X(L.J) 
X(L.J)=X(KP.J) 

15  X(KP.J)=Z 

20  IF{ABSF( A(L.L) )-EP)50. 50.30 

30  IF(L-N)31 ,34,34 

31  LPl=l+l 

DO  36  K=LP1,N 
IF(A(K,L) )32,36,32 

32  RATIO=A(K,L)/A(L,L) 
DO  33  J=LP1,N 

33  A(K,J)=A(K.J)-RATIO*A(L»J) 
DO  35  J=1,N 

35  X(K»J)=X(K,J)-RATIO*X(L.J) 

36  CONTINUE 

34  CONTINUE 

40  DO  43  I=1.N 
II=N+1-I 

DO  43  J=l .N 

S  =  0.0 

IF( I I-N)41,43»A3 

41  I  I P 1  =  I  I  +  1 

DO  42  K=I IP1»N 

43  X{ I I»J)  =  (X( I  I »J)-S)/A( II. I  I ) 
KER  =  1 
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RETURN 
50  KER=2 

END 

SUBROUTINE  PHIDEL  ( PH I , DEL , N , DT ) 

DIMENSION  F(12, 12), PH 1(12,12) ,TFRM(12,12) ,WORM( 12,12) 
1  ,  DEL(  12)  ,DELM(  12*12)  »TELMU2,12)  »DELP(  12»12),D(12) 
C      VALID  ONLY  FOR  A  CONSTANT  F  MATRIX 
C      DT=  SAMPLING  INTERVAL 

NFINAL=1 

READ  1(  (F(  IR,  IC)  ,IC=1,N) ,IR=1,N> 

READ  1  <D( I  )  ,1  =  1, N) 

1  FORMAT  ((8F10.0)) 
TM=0.0 

1003  PRINT  399, DT,  (  ( F( IR,  IC)  ,IC  =  1,N)  ,IR  =  1,N) 

399  F0RMAT(///4H  DT=,  1F5.3////,8H  F ( I , J ) =/ , 4 ( 4E 16. 7/ ) ) 
PRINT  3991  (D(  I  )  ,  I  =  1  ,  IM  ) 

3991  FORMAT(///6H  D (  I ) = / , 4 ( 1 E16  .7/ )  ) 
DO  400  IR=1,N 
DO  400  IC=1,N 
TERM( I R  » IC ) =0.0 
WORM( IR,IC) =0.0 
TERM( IR,IR) =1.0 
TELM( IR,IC)=TERM( IR, IC)*DT 
DELP( IR,IC) =TELM{ IR, IC) 
DELM( IR,IC) =0.0 
DEL( IR)=0.0 

400  PHI ( IR, IC)=TERM( IR,IC) 
4  TM=1.0+TM 

DO  5  00  IR=1,N 
DO  500  IC=1,N 
DO  500    JN=1,N 

DELM( IR, IC) =DELM( IR, IC)-TELM( IR,JN)*F( JN,IC)*DT/(TM+1.0) 
500  WORM( IR, IC) =TERM( I R, JN ) *F ( JN , I C ) *DT/TM+WORM( IR,IC) 
DO  401  IR=1,N 
DO  401  IC=1,N 
TERM( IR,IC) =WORM( IR, IC  ) 
TELM( I R , IC ) =DELM( IR, IC) 
DELP(  IR,  IC) =DELP(  I R ,  I C ) +TELM ( I R  ,  I C ) 
PHI(IR,IC)=PHI(IR,IC)+TERM(IR,IC) 
DFLM( IR,IC) =0.0 

401  WORM( IR,IC) =0.0 
ABC=0.0 
DO  21  =  1, N 
DO  2J=1,N 
AA=TERM( I ,J) 
AB=ABSF( AA) 
IF(ABC-AB)3,3,2 

3  ABC=A3 

2  CONTINUE 
IF(0.0  000  00n0  5-ARC ) 5  ,5,6 

a  ^fi  tq  ft 

6  PRINT  502, (  (PHI  (  IR,IC)  ,  IC  =  1,N)  , IR  =  1,N) 
50  2  FORMAT ( //9X,8HPHI(I,J)///((4El6.7))) 
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8  £83 

DO  600 


J=1,N 

60  0  DFL(  I  )=DEL(  I  )+PHI(  I  ,  J)  *DELP  (  J  ,K  )  *D<  K  » 

PRINT  503  (DELC I  )  ,  I  =  1»N> 
50  3  FORMAT {///9X»6HDEL( I)//(8F15.Q)//) 

RETURN 

END 

FND 

• 
PROGRAM  MISSILE  CONTINUED 


EVASION 


YTDOT  =YTDOT+10. 

EVASION  LEFT-RIGHT-LEFT 
IF(<K-3O)320, 321,321 

320  YTDOT=YTDOT+20. 
GOTO  323 

321  YTDOT  =  YTDOT-20.  . 

323  CONTINUE 
IF(YTZ-5O00. ) 324,324,32  5 

324  YTDOT=YTDOT+30. 

IF  <YTZ-2000.)326,326,325 
326  YTDOT=0.0 

325  CONTINUE 

END  OF  EVASION 


a7 


APPENDIX  II 

The  missile  guidance  simulation  is  a  transfer  function  relating 
(  y  )  and  (  a  ) .    Since  the  objective  of  this  study  is  to  examine  the 
effect  of  a  digital  filter,  an  arbitrary  transfer  function  was  se- 
lected which  could  represent  any  present  day  missile.    This 
transfer  function  contains  a  time  delay  due  to  the  radar,  a  fil- 
tering action  due  to  the  control  (hydraulics  etc.),  and  a  second 
order  system  representing  the  missile  dynamics. 

°/y  =         A  .  B  .  C  (A-l) 

S  +  TR  S  +  T  Sz  +  DS  +  E 

As  mentioned  above  the  intention  here  is  to  demonstrate  the 
capability  of  the  digital  filter  in  the  presence  of  a  large  mag- 
nitude of  noise.    The  values  of  the  gains  and  time  constants  were 
selected  for  equation  (A-l)  after  examining  several  different 
missile  systems  which  use  proportional  navigation.    Refer  to 
[13]  for  further  elaboration  on  determining  the  F  and  D  matrices 
for  a  particular  system . 
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